Data prep

Cleaning and formatting

We start by preparing the data from the outputs of scripts 7.1 and 8.1. The results table shown as output is the key to read the figures.

# Results dir + mun shapefile
mun <- st_read("../data/mun/munic_SHP_clean.shp", quiet = TRUE)

# Classes
classes <- read_csv("../config/rcl_tables/land_use/recode_table.csv")
forest_classes <- read_csv("../config/rcl_tables/land_use/recode_table_forest.csv") %>% 
  rename(new_code = "ID", new_class = "Name")

classes_unique <- unique(classes[, c("new_code", "new_class")]) %>% 
  bind_rows(forest_classes)

# results key
sce_code_vec <- 
  as.vector(t(outer(c("BAU-", "R-", "CorPr-", "R-CorPr-", "R(T)-CorPr-"), 
                    c("Hist", "Base", "RCP8"), paste0)))
sce_code_vec_run <- 
  rep(c("BAU", "R", "CorPr", "R-CorPr", "R(T)-CorPr"), each=3)

results_clean <- read_rds("../data/temp/stsim_run_results.RDS") %>% 
  dplyr::select(scenarioId, name) %>% 
  mutate(name =  gsub(.$name, pattern = " \\(.+\\)", replacement = "")) %>% 
  mutate(sce =  paste0("sce_", .$scenarioId)) %>% 
  mutate(chapter = c("none", rep("both", 3), rep("chap_1", 3), rep("chap_2", 9))) %>% 
  mutate(splitted = str_split(name, " \\| ")) %>% 
  mutate(climate = unlist(map(splitted, ~unlist(.x[2])))) %>% 
  mutate(run = unlist(map(splitted, ~unlist(.x[1])))) %>% 
  dplyr::select(-splitted) %>% 
  replace_na(list(climate = "none")) %>% 
  mutate(code = c("Control", sce_code_vec)) %>% 
  mutate(code_run = c("control", sce_code_vec_run))

# Final extracted datasets
df_final <- readRDS("../outputs/final/final_df_current_density.RDS") %>%
  mutate(timestep = (timestep*10)+1990, source = "model")
df_final_origin <- readRDS("../outputs/final/final_df_origin_current_density.RDS") %>%
  mutate(timestep = timestep*10+1990, source = "model")

# Stratum key
stratum_key <- read_csv("../config/stsim/SecondaryStratum.csv") %>%
  mutate(ID = as.factor(ID)) %>%
  rename(zone=ID, MUS_NM_MUN=Name) %>%
  left_join(df_final, by="zone") %>%
  filter(MUS_NM_MUN!="Not Monteregie")

# Sum all municipalities
df_summarised <- df_final %>%
  group_by(sce, timestep, species, iteration) %>% 
  summarise(sum_cur = sum(current)) %>% ungroup %>%
  mutate(source = "model")
df_origin_summarised <- df_final_origin %>%
  group_by(timestep, species) %>% 
  summarise(sum_cur = sum(mean)) %>% ungroup %>% 
  mutate(sce = "sce_0", source = "observation")

# Joined dataset
joined <- full_join(df_summarised, df_origin_summarised, 
                    by=c("sce", "source", "species", "timestep", "sum_cur")) %>% 
  left_join(results_clean, by = "sce")  %>% 
  replace_na(list(climate = "none", run = "historic run")) %>% 
  # Make factors
  mutate(sce = as.factor(sce), run = as.factor(run), sum_cur = 10*(sum_cur),
         climate = factor(climate, levels = c("none", "historic",
                                              "baseline", "RCP 8.5")),
         run = factor(run, levels = c("historic run", "BAU run", 
                                      "BAU run + corrs protection", "BAU run + ref", 
                                      "BAU run + corrs protection + ref",
                                      "BAU run + corrs protection + ref TARGETED")))

# Roc curves data
urb <- readRDS("../data/temp/fit_rs_outcome_rf_urb_2.RDS")
agex <- readRDS("../data/temp/fit_rs_outcome_rf_agex_2.RDS")

# Histogram data 
hist_original <- read_csv("../outputs/final/final_values_output_original.csv")
hist_true <- read_csv("../outputs/final/final_values_output_TRUE.csv") %>% 
  mutate(sce = "TRUE")
histograms <- read_csv("../outputs/final/final_values_output.csv") %>% 
  bind_rows(hist_original) %>% 
  bind_rows(hist_true) %>% 
  left_join(results_clean, by="sce") 

# SURF Data
surf <- read_csv("../surf/surf_output.csv") %>% 
  mutate(timestep = timestep*10+1990) %>% 
  left_join(results_clean, by=c("scenario"="scenarioId")) %>% 
  replace_na(list(climate = "none", run = "historic run")) %>% 
  # Make factors
  mutate(sce = as.factor(sce), run = as.factor(run),
         climate = factor(climate, levels = c("none", "historic",
                                              "baseline", "RCP 8.5")),
         run = factor(run, levels = c("historic run", "BAU run", 
                                      "BAU run + corrs protection", "BAU run + ref", 
                                      "BAU run + corrs protection + ref",
                                      "BAU run + corrs protection + ref TARGETED")))

# Bar plot data
bar_data <- read_csv("../outputs/final/final_bar_plot_data.csv")

results_clean

Figures

Mean current flow change

Historic

the_width = 18 
the_height = 15
the_dpi = 500

fig_1_historic <- joined %>% 
  filter(climate == "none") %>% 
  group_by(scenarioId, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 1990)$sum_cur, 
                       the_diff = ((sum_cur-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=source)) +
  scale_color_manual(values=c('#d8b365','#5ab4ac'), 
                     labels = c("Model", "Observation")) +
  geom_smooth(aes(group = sce), method = "glm", se=FALSE) +
  scale_linetype_manual(values = c(1,3,2,5))+
  geom_point(aes(group = sce)) +
  facet_wrap(~species, scales = "fixed") +
  labs(y = "Current flow (% Change)",
       x = "Year",
       col = "Source") +
  NULL
fig_1_historic
Connectivity change for species through time, 1990-2010

Connectivity change for species through time, 1990-2010

ggsave(fig_1_historic,
       filename = "../thesis/figures/connectivity_decrease_x5species_historic.png",
       width = the_width, height = the_height, dpi = the_dpi)

Chap 1

Line graph
colors_sce <- 
  data.frame(sce = c("Hist", "Baseline", "RCP85"), 
             color = c("#8DA0CB", "#FC4E2A", "#800026"), stringsAsFactors = F)

fig_1_static_chap_1 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(scenarioId, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$sum_cur, 
                       the_diff = ((sum_cur-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color, 
                     labels = c("Historic", "Baseline", "RCP 8.5")) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  #geom_point(aes(group = sce)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  #geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "Current flow (% Change)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL
fig_1_static_chap_1
Connectivity change for species through time, 2010-2100

Connectivity change for species through time, 2010-2100

ggsave(fig_1_static_chap_1, 
       filename = "../thesis/figures/connectivity_decrease_x5species_chap1.png", 
       width = the_width, height = the_height, dpi = the_dpi)
Radar Graph
# Default size
gridline.label.offset   = 4
legend.text.size = 20
group.line.width = 0.9
grid.label.size = 10
background.circle.colour    = "#ffffff"
group.point.size    = 5
axis.label.size = 8
group.colours = RColorBrewer::brewer.pal(name = "Paired",n=10)[c(2,4,6,8,10)]


radar_data_chap1 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  group_by(timestep, species, code) %>% 
  summarise(sum_cur = mean(sum_cur)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = sum_cur) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

radar_chap1 <-  
  ggradar(radar_data_chap1, centre.y = -20, legend.position = "right",
          grid.min = -20, grid.max = 3, grid.mid = 0, 
          values.radar = c("-20 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = gridline.label.offset, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size)
radar_chap1

ggsave("../thesis/figures/radar_ggradar_chap1.png", radar_chap1,
       width = the_width, height = the_height, dpi = the_dpi)

Chap 2

fig_1_static_chap_2 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  group_by(scenarioId, species, iteration) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$sum_cur, 
                       the_diff = ((sum_cur-baseline)/baseline)*100 )) %>% 
  ggplot(aes(x=timestep, y=the_diff, col=climate)) +
  scale_color_manual(values = colors_sce$color,
                     labels = c("Historic", "Baseline", "RCP 8.5")) +
  geom_smooth(aes(group = sce, linetype=code_run), alpha=0.5, method="loess") +
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed") +
  geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 2020, linetype = "dashed", color = "darkGreen", alpha = 0.5) +
  labs(y = "Current flow (% Change)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow = 2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  theme(legend.position = c(0.85, 0.15))+
  NULL
fig_1_static_chap_2
Connectivity change for species through time, 2010-2100

Connectivity change for species through time, 2010-2100

ggsave(fig_1_static_chap_2, 
       filename = "../thesis/figures/connectivity_decrease_x5species_chap2.png", 
       width = the_width, height = the_height, dpi = the_dpi)
Radar Graph
radar_data_chap2 <- joined %>% 
  filter(climate != "none", chapter %in% c("none", "both","chap_2")) %>% 
  group_by(timestep, species, code) %>% 
  summarise(sum_cur = mean(sum_cur)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = sum_cur) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

radar_chap2 <-  
  ggradar(radar_data_chap2, centre.y = -20, legend.position = "right",
          grid.min = -20, grid.max = 3, grid.mid = 0, 
          values.radar = c("-20 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = gridline.label.offset, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size-2)
radar_chap2

ggsave("../thesis/figures/radar_ggradar_chap2.png", radar_chap2,
       width = the_width, height = the_height, dpi = the_dpi)

Both

Radar Graph
radar_data_all <- joined %>% 
  filter(climate != "none") %>% 
  group_by(timestep, species, code) %>% 
  summarise(sum_cur = mean(sum_cur)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = sum_cur) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

radar_all <-  
  ggradar(radar_data_all, centre.y = -20, legend.position = "right",
          grid.min = -20, grid.max = 3, grid.mid = 0, 
          values.radar = c("-20 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = gridline.label.offset, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size-2)
radar_all

ggsave("../thesis/figures/radar_ggradar_both.png", radar_all,
       width = the_width, height = the_height, dpi = the_dpi)

ROC Curves

ROC curves from fitting both models, plotting with patchwork package.

p1 <- bind_rows(urb, .id = "fold") %>% 
  mutate(Fold = factor(fold, levels = as.character(1:10))) %>% 
  group_by(Fold) %>% 
  roc_curve(.pred, truth = outcome_fact) %>% 
  autoplot(add=T) +
  ggtitle(paste("Urbanisation")) +
  annotate(x = 0.75, y = 0.25, geom="label", size = 3,
           label = as.character(paste("Av AUC = 0.938 +/- 0.002"))) +
  theme(legend.position = "none") +
  NULL

p2 <- bind_rows(agex, .id = "fold") %>% 
  mutate(Fold = factor(fold, levels = as.character(1:10))) %>% 
  group_by(Fold) %>% 
  roc_curve(.pred, truth = outcome_fact) %>% 
  autoplot(add=T) +
  ggtitle(paste("Agricultural Expansion")) +
  annotate(x = 0.75, y = 0.25, geom="label", size = 3,
           label = as.character(paste("Av AUC = 0.929 +/- 0.002"))) +
  NULL

full = p1 + p2

full

ggsave("../thesis/figures/double_roc_resample.png", full)

Histograms

Original

bi_colors <- c("#67a9cf", "#ef8a62")

hist_plot_origin <- histograms %>%
  filter(sce %in% c("TRUE", "sce_37")) %>% 
  filter(ts == 2) %>% 
  mutate(ts = ts*10+1990) %>% 
  mutate(ts = as.factor(ts)) %>% 
  ggplot(aes(x=bins, y=n, fill=sce, colour =sce)) + 
  geom_density(stat='identity', show.legend=T, alpha=0.3)+
  scale_fill_manual(values = bi_colors, name="Source", labels=c("Model", "Observation")) +
  scale_color_manual(values = bi_colors, name="Source", labels=c("Model", "Observation")) +
  facet_wrap(~species) +
  labs(x = "Flow intensity distribution (log)", 
       y = "") +
  theme(strip.text.y = element_text(angle=360, size=10, hjust = 0), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.text.y.left = element_blank(),
        axis.ticks.y.left = element_blank()) +
  NULL
hist_plot_origin

ggsave("../thesis/figures/original_hists.png", hist_plot_origin,
       width = the_width, height = the_height, dpi = the_dpi)

Chap 1

bi_colors <- c("#67a9cf", "#ef8a62")

hist_plot_1 <- histograms %>% 
  mutate(ts = as.factor(ts)) %>% 
  filter(chapter %in% c('both', 'chap_1')) %>%  
  ggplot(aes(x=bins, y=n, group=ts)) + 
  geom_density(stat='identity', show.legend=T, aes(fill=factor(ts), color=factor(ts)), alpha=0.3)+
  facet_grid(code~species) +
  scale_fill_manual(values = bi_colors, 
                    labels = c("2010", "2100")) +
  scale_color_manual(values = bi_colors, 
                     labels = c("2010", "2100")) +
  labs(x = "Flow intensity distribution (log)", 
       y = "")+
  theme(strip.text.y = element_text(angle=360, size=10, hjust = 0), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.text.y.left = element_blank(),
        axis.ticks.y.left = element_blank(), 
        legend.position = c(0.95, 0.96)) +
  labs(color="Timestep", fill="Timestep") +
  NULL
hist_plot_1

ggsave("../thesis/figures/hist_chap1.png", hist_plot_1, width = 18, height = 15)

Chap 2

hist_plot_2 <- histograms %>% 
  mutate(ts = as.factor(ts)) %>% 
  filter(chapter %in% c('chap_2')) %>%  
  ggplot(aes(x=bins, y=n, group=ts)) + 
  geom_density(stat='identity', aes(fill=factor(ts),  colour=factor(ts)), alpha=0.3)+
  facet_grid(code~species) +
  scale_fill_manual(values = bi_colors, 
                    labels = c("2010", "2100")) +
  scale_color_manual(values = bi_colors, 
                     labels = c("2010", "2100")) +
  labs(x = "Flow intensity distribution (log)", 
       y = "") +
  theme(strip.text.y = element_text(angle=360, size=10, hjust = 0), 
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(), 
        axis.text.y.left = element_blank(),
        axis.ticks.y.left = element_blank(),
        legend.position = c(0.95, 0.96)) +
  labs(color="Timestep", fill="Timestep") +
  NULL
hist_plot_2

ggsave("../thesis/figures/hist_chap2.png", hist_plot_2, width = 18, height = 20)

SURF Analysis

Chap1

Linear Graph
surf_1 <- surf %>% 
  mutate(scenario = as.factor(scenario)) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_1")) %>% 
  filter(climate != "none") %>% 
  group_by(scenario, species, iter) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$kp_nb, 
                       the_diff = ((kp_nb-baseline)/baseline)*100 )) %>% ungroup %>% 
  ggplot(aes(x = timestep, y = the_diff, color=climate)) +
  scale_linetype_manual(values=c("solid", "twodash")) +
  geom_smooth(aes(group = scenario, linetype=code_run), alpha=0.2, method="loess")+
  #geom_point(aes(group=scenario), alpha= 0.05, size = 0.5)+
  facet_wrap(~species, scales = "fixed")+
  scale_color_manual(values = colors_sce$color,
                     labels = c("Historic", "Baseline", "RCP 8.5")) +
  labs(y = "Change in nb of Keypoints detected (%)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=1, override.aes = list(colour = 'black'))) +
  NULL
surf_1

ggsave("../thesis/figures/surf_chap1.png", surf_1,
       width = the_width, height = the_height, dpi = the_dpi)
Radar Graph
surf_radar_data_1 <- surf %>% 
  filter(climate != "none", chapter %in% c("none", "both","chap_1")) %>% 
  filter(sce != "sce_0", name != "historic run") %>% 
  group_by(timestep, species, code) %>% 
  summarise(kp_nb = mean(kp_nb)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = kp_nb) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

surf_radar_1 <- 
  ggradar(surf_radar_data_1, centre.y = -60, legend.position = "right",
          grid.min = -60, grid.max = 5, grid.mid = 0, 
          values.radar = c("-60 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = 20, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size)
surf_radar_1

ggsave("../thesis/figures/surf_radar_chap1.png", surf_radar_1,
       width = the_width, height = the_height, dpi = the_dpi)

Chap 2

Linear Graph
surf_2 <- surf %>% 
  mutate(scenario = as.factor(scenario)) %>% 
  filter(climate != "none", chapter %in% c("none", "both", "chap_2")) %>% 
  filter(climate != "none") %>% 
  group_by(scenario, species, iter) %>% 
  group_modify(~mutate(.x , baseline = subset(.x, timestep == 2010)$kp_nb, 
                       the_diff = ((kp_nb-baseline)/baseline)*100 )) %>% ungroup %>% 
  ggplot(aes(x = timestep, y = the_diff, color=climate)) +
  geom_smooth(aes(group = scenario, linetype=code_run), alpha=0.2, method="loess")+
  #geom_point(aes(group=scenario), alpha= 0.05, size = 0.5)+
  scale_linetype_manual(values=c("solid", "dashed", "longdash", "dotted")) +
  facet_wrap(~species, scales = "fixed")+
  scale_color_manual(values = colors_sce$color,
                     labels = c("Historic", "Baseline", "RCP 8.5")) +
  labs(y = "Change in nb of Keypoints detected (%)",
       x = "Year",
       col = "Climate",
       linetype = "Run") +
  guides(col = guide_legend(ncol = 2, nrow=2), 
         linetype = guide_legend(nrow=2, ncol=2, 
                                 override.aes = list(colour = 'black'))) +
  NULL
surf_2

ggsave("../thesis/figures/surf_chap2.png", surf_2,
       width = the_width, height = the_height, dpi = the_dpi)
Radar Graph
surf_radar_data_2 <- surf %>% 
  filter(climate != "none", chapter %in% c("none", "both","chap_2")) %>% 
  filter(sce != "sce_0", name != "historic run") %>% 
  group_by(timestep, species, code) %>% 
  summarise(kp_nb = mean(kp_nb)) %>% ungroup %>% 
  pivot_wider(names_from = timestep, values_from = kp_nb) %>% 
  mutate(diff = ((`2100`-`2010`)/`2010`)*100) %>% 
  select(-c(`2010`:`2100`)) %>%  
  pivot_wider(names_from = code, values_from = diff) %>% 
  rename(group = species)

surf_radar_2 <- 
  ggradar(surf_radar_data_2, centre.y = -60, legend.position = "right",
          grid.min = -60, grid.max = 5, grid.mid = 0, 
          values.radar = c("-60 %", "0 %", ""), 
          group.colours = group.colours,
          gridline.label.offset = 22, 
          legend.text.size = legend.text.size,
          group.line.width = group.line.width, 
          grid.label.size = grid.label.size,
          background.circle.colour = background.circle.colour, 
          group.point.size  = group.point.size, 
          axis.label.size = axis.label.size-2)
surf_radar_2

ggsave("../thesis/figures/surf_radar_chap2.png", surf_radar_2, 
       width = the_width, height = the_height, dpi = the_dpi)

Maps

Land use maps

source('~/Documents/Master/Thesis/land_con_monteregie/scripts/functions/draw_scenario.R')
for (scenario in c(39, 42, 45, 48, 51)){
  the_plot <- draw_scenario(scenario, mxp = mxp)
  the_name <- paste0("scenario_", as.character(scenario),"_compare.png")
  ggsave(plot = the_plot, 
         filename = paste0("../thesis/figures/",the_name), width = 20, height = 10, dpi = the_dpi)
}
## Parsed with column specification:
## cols(
##   source = col_character(),
##   name = col_character(),
##   class = col_character(),
##   subclass = col_logical(),
##   old_code = col_double(),
##   new_code = col_double(),
##   new_class = col_character()
## )
## Parsed with column specification:
## cols(
##   ID = col_double(),
##   Name = col_character()
## )
## Parsed with column specification:
## cols(
##   source = col_character(),
##   name = col_character(),
##   class = col_character(),
##   subclass = col_logical(),
##   old_code = col_double(),
##   new_code = col_double(),
##   new_class = col_character()
## )
## Parsed with column specification:
## cols(
##   ID = col_double(),
##   Name = col_character()
## )
## Parsed with column specification:
## cols(
##   source = col_character(),
##   name = col_character(),
##   class = col_character(),
##   subclass = col_logical(),
##   old_code = col_double(),
##   new_code = col_double(),
##   new_class = col_character()
## )
## Parsed with column specification:
## cols(
##   ID = col_double(),
##   Name = col_character()
## )
## Parsed with column specification:
## cols(
##   source = col_character(),
##   name = col_character(),
##   class = col_character(),
##   subclass = col_logical(),
##   old_code = col_double(),
##   new_code = col_double(),
##   new_class = col_character()
## )
## Parsed with column specification:
## cols(
##   ID = col_double(),
##   Name = col_character()
## )
## Parsed with column specification:
## cols(
##   source = col_character(),
##   name = col_character(),
##   class = col_character(),
##   subclass = col_logical(),
##   old_code = col_double(),
##   new_code = col_double(),
##   new_class = col_character()
## )
## Parsed with column specification:
## cols(
##   ID = col_double(),
##   Name = col_character()
## )

Bar charts

Sce 39

final_df_39 <- bar_data %>% 
  filter(scenario == 39)

bar_plot_39 <- final_df_39 %>% 
  #filter(value > 10) %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  #geom_errorbar(aes(ymin=mean_count-sd_count, ymax=mean_count+sd_count), width=.2,
  #               position=position_dodge(.9)) +
  scale_fill_manual(values = bi_colors) +
  theme(legend.position = c(0.85, 0.7)) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL
bar_plot_39

#ggsave("../thesis/figures/bar_BAU.png", bar_plot_39, width = 15, height = 10)
Sce 42
# draw_scenario(42)
final_df_42 <- bar_data %>% 
  filter(scenario == 42)

bar_plot_42 <- final_df_42  %>% 
  #filter(value > 10) %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(values = bi_colors) +
  #geom_errorbar(aes(ymin=mean_count-sd_count, ymax=mean_count+sd_count), width=.2,
  #               position=position_dodge(.9)) +
  theme(legend.position = c(0.85, 0.7)) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL
bar_plot_42

#ggsave("../thesis/figures/bar_Ref.png", bar_plot_42, width = 15, height = 10)
Both
final_df_39 <- final_df_39 %>% 
  mutate(sce = "BAU") %>% 
  mutate(type = ifelse(value <10, "non_forest", "forest"))
final_df_42 <- final_df_42 %>% 
  mutate(sce = "Reforestation")  %>% 
  mutate(type = ifelse(value <10, "non_forest", "forest"))

final_all <- bind_rows(final_df_39, final_df_42)

forest <- final_all %>% 
  filter(type == "forest") %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(values = bi_colors) +
  theme(legend.position = "none") +
  facet_grid(~sce) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL

non_forest <- final_all %>% 
  filter(type != "forest") %>% 
  mutate(value = as.factor(value), timestep = as.factor(timestep)) %>% 
  ggplot(aes(fill=timestep, y=area, x=new_class)) + 
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(values = bi_colors) +
  theme(legend.position = "right") +
  facet_grid(~sce) +
  labs(x = "Land use Class", 
       y = "Mean Area (hectare)", 
       fill = "Year") +
  NULL

assembled <- wrap_plots(forest, non_forest, ncol=1, nrow=2)
assembled

ggsave("../thesis/figures/bar_both.png", assembled,
       width = the_width, height = 20, dpi = the_dpi)